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Dielectric optical micro-resonators and micro-lasers represent a realization of a wave-chaotic sys- 
tem, where the lack of symmetry in the resonator shape leads to non-integrable ray dynamics. 
Modes of such resonators display a rich spatial structure, and cannot be classified through mode 
indices which would require additional constants of motion in the ray dynamics. Understanding and 
controlling the emission properties of such resonators requires the investigation of the correspon- 
dence between classical phase space structures of the ray motion inside the resonator and resonant 
solutions of the wave equations. We first discuss the breakdown of the conventional eikonal approx- 
imation in the short wavelength limit, and motivate the use of phase-space ray tracing and phase 
space distributions. Next, we introduce an efficient numerical method to calculate the quasi-bound 
modes of dielectric resonators, which requires only two diagonalizations per N states, where N is 
approximately equal to the number of half-wavelengths along the perimeter. The relationship be- 
tween classical phase space structures and modes is displayed via the Husimi projection technique. 
Observables related to the emission pattern of the resonator are calculated with high efficiency. 

I. INTRODUCTION 

A promising approach to making compact and liigh-Q optical resonators is to base them on "totally internally 
reflected" modes of dielectric micro-structures. Such devices have received considerable attention as versatile com- 
ponents for integrated optics and for low threshold micron-scale semiconductor lasers 0, l6^ . The interest in such 
resonators for applications and for fundamental optical physics has motivated the extension of optical resonator 
theory to describe such systems. 

All optical resonators are open systems described by modes characterized by both a central frequency and a width 
(their ratio giving the Q-factor of the mode) . In a mirror-based resonator the set of resonant frequencies is determined 
by an optical path-length for one round-trip along a path determined by the mirrors within the resonator; the width 
is determined by the reflectivity of the mirrors, diffraction at the mirror edges and by absorption loss within the 
resonator. Accurate analytic formulas can be found for the resonator frequencies and for the electric field distribution 
of each mode using the methods of Gaussian optics [g^- The modes are characterized by one longitudinal and 
two transverse mode indices (in three dimensions). These mode indices play the same role mathematically for the 
electromagnetic wave equation as the good quantum numbers play in characterizing solutions of the wave equation of 
quantum mechanics. 

For fabricating optical resonators on the micron scale, using total internal reflection from a dielectric interface for 
optical conflnement is convenient as it simplifies the process. Such dielectric resonators deflne no specific optical path 
length; many different and potentially non-closed ray trajectories can be confined within the resonator. An important 
point, emphasized in the current work, is that in such resonators there typically exist many narrow resonances 
characterized by their frequency and width, but such resonances often cannot be characterized by any further modes 
indices. This is the analog of a quantum system in which there are no good quantum numbers except for the energy. 
We shall see that the way to determine if a given mode has additional mode indices (other than the frequency), is 
to determine whether it corresponds to regular or chaotic ray motion. We will present below an efficient numerical 
method for calculating all of the resonances of a large class of dielectric resonators; we will also describe the surface 
of section and Husimi- Wigner projection method to determine the ray dynamics corresponding to such a mode. 

Although both DBR-based and edge-emitting optical resonators rely on reflectivity from a dielectric interface (at 
normal incidence), we will use the term dielectric resonator (DR) to refer to resonators that rely on the high reflectivity 
of dielectric bodies to radiation incident from within the dielectric near the critical angle for total internal reflection. 
This is the only class of resonators we will treat below. We immediately point out that totally-internally-reflected 
solutions of the wave equation only exist for inflnite flat dielectric interfaces; any curvature or finite extent of the 
dielectric will allow evanescent leakage of propagating radiation from the optically more dense to the less dense 
medium. As a dielectric resonator is a finite dielectric body embedded in air (or in a lower index medium) it will of 
necessity allow some evanescent leakage of all modes, even those which from ray analysis appear to be totally-internally 
reflected. 

A very large range of shapes for DRs have been studied during the recent years. By far the most widely studied are 
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rotationally symmetric structures such as spheres, cyhnders and disks. The reason for this is that the wave equation 
is separable and the solutions can be written in terms of special functions carrying three modes indices. The narrow 
(long-lived) resonances correspond to ray trajectories circling around the symmetry axis near the boundary with angle 
of incidence above total internal reflection; these solutions are often referred to as "whispering gallery" (WG) modes 
or morphology-dependent resonances. In this case, due to the separability of the problem, it is straightforward to 
evaluate the violation of total internal reflection, which may be interpreted as the tunneling of waves through the 
angular momentum barrie r l33l l48j . Micron-scale, high-Q micro-lasers were fabricated in mid-80s and early 90s based 
on such cylindrical ^3iE3E3 (disk-shaped) {Q ^ 10* — 10^) and spherical [l^ {Q ~ 10^ — 10^^) dielectric resonators. 
However, the very high Q value makes these resonators unsuitable for micro-laser applications, because such lasers 
invariably provide low-output power and furthermore, unless additional guiding elements are used, the lasing output 
is emitted isotropically. 

As early as 1994|53 one of the current authors proposed to study dielectric resonators based on smooth deformations 
of cylinders or spheres which were referred to as "asymmetric resonant cavities" (ARCs) . The idea was to attempt 
to combine the high Q provided by near total internal reflection with a breaking of rotational symmetry leading to 
directional emission and improved output coupling. General principles of non-linear dynamics applied to the ray 
motion (to be reviewed below) suggested that there would be only a gradual degradation of the high-Q modes, and 
one might be able to obtain directional emission from deformed whispering gallery modes. Expcrimental |lOl l46l l5l| 
and theoretical^^ work since that initial suggestion has confirmed this idea, although the important modes are not 
always of the whispering gallery tvpe,21, 22, 23, 40, 59j. 

Calculating the modal properties of deformed cylindrical and spherical resonators presents a much more challenging 
theoretical problem. Unless the boundary of the resonator corresponds to a constant coordinate surface of some or- 
thogonal coordinate system, the resulting partial differential equation will not be solvable by separation of variables. 
The only relevant separable case is an exactly elliptical deformation of the boundary, which turns out to be unrepre- 
sentative of generic smooth deformations. Using perturbation theory to evaluate the new modes based on those of the 
cylindrical or spherical case is also impractical, as for interesting deformations and typical resonator dimensions (tens 
of microns or larger) the effect of the deformation is too large for the modes of interest to be treated by perturbation 
theory. The small parameter in the problem for attempting approximate solutions is the ratio of the wavelength 
to the perimeter X/2t:R = (kR)^^ . Eikonal methods {Sq and Gaussian optical methods |63| both rely on the short 
wavelength limit to find approximate solutions. The Gaussian optical method can be used to find a subset of the 
solutions for generic ARCs, those associated with stable periodic ray orbits, as explained in detail in Ref . [65l| . The 
eikonal method can also be used to find a subset of the modes of ARCs, if one has a good approximate expression for 
a local constant of motion; an example of this is the adiabatic approximation used by Nockel and Stone ji^. However 
a large fraction of the modes in ARCs are not describable by either of these methods. The breakdown of the Gaussian 
optical methods is easily seen as a fully chaotic system will have only unstable periodic orbits and the solutions one 
obtains by the Gaussian method near unstable periodic ray orbits are inconsistent |65l| . The failure of eikonal methods 
is more subtle and really arises from the possibility of chaotic ray motion in a finite fraction of the phase space. 
Often optics textbooks and even standard research references treat the eikonal method as being of completely general 
applicability; we therefore will devote the next section of this paper to an explanation of the failure of eikonal methods 
for resonators with arbitrary smooth boundaries. In section (|III|1 we describe the phase space methods which indicate 
that this failure is generic. 

In section (|IV|) we present the formulation of the resonance problem and in section the reduction of the Maxwell's 
equations to the Helmholtz Equation for the resonators we study. The failure of all standard short wavelength 
approximation methods to describe the solution of wave equations in finite domains with arbitrary smooth boundaries 
has led to the problem of "quantizing chaos" in the context of the Schrodingcr equation and the Helmholtz Equation 
(although our problem is somewhat different due to the dielectric boundary conditions on this equation). Although 
substantial progress has been made using periodic orbit methods in obtaining approximations for the density of states 
of fully chaotic systems, these methods do not yield individual solutions of the wave equation. It is therefore of 
great importance in this field to develop efhcient numerical methods which can be used to calculate and interpret 
the resonance properties. In section ljVI|l we present a highly efficient numerical method for ARCs, adapted from the 
S-matrix methods developed in the field of quantum chaos. In section (|VIII| we display a range of resonant solutions 
for a partially chaotic dielectric resonator and in section ljVIII|) we show how to perform the Husimi projection of the 
real-space numerical solutions so obtained into phase space in order to interpret them in terms of ray dynamics. The 
calculation of experimental observable relating to emission patterns from micro-lasers are discussed in section (|IXI| . 
Finally, in section lpC)l . we show examples of the main types of modes one encounters in wave-chaotic dielectric 
resonators generally and specifically in ARCs. 
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II. FAILURE OF EIKONAL METHODS FOR GENERIC DIELECTRIC RESONATORS 

The use of classical ray theory to describe monochromatic, high-frequency solutions of the waye equation is described 
in yarious references [3. l35l l38|| . The connection between rays and waycs is standardly deriycd in the context of the 
Helmholtz equation 

{V^ + n^{x)k^)^{x) ^0; (1) 

the waye equation for the resonator problem will be reduced to this equation in section l(V)l . The eikonal approach 
uses the asymptotic ansatz 

^(.) ^ e'^^-^t^) Y: ^ (2) 

1^=0 

in the limit k — > oo. Inserting Eq. Q into Eq. one finds to lowest order in the asymptotic parameter 1/k, the 
eikonal equation 

[VSf = n\x) (3) 

and the transport equation 

2VS -VAq + AoW^S ^0 (4) 

Note at this point we only assume one eikonal S{x) and one amplitude A{x) at each order in the expansion. We will 
also specialize to a uniform medium of dielectric constant n. In this framework, each waye solution 'ip{x) corresponds 
to a family of rays defined by the ycctor field 

p{x) = VS{x) (5) 

where the field has a fixed magnitude, IVS*] = n. The solution for the function S{x) can be found by the specification 
of initial value boundary conditions on an open curye C : x = x{s) and propagating the curve using the eikonal 
equation. Such an initial value solution can thus be extended until it encounters a point at which two or more distinct 
rays of the wavefront converge; at or nearby such a point will occur a focus or caustic at which the amplitude A will 
diverge and in the neighborhood of which the asymptotic representation becomes ill-defined. (A c austic is a curve to 
which all the rays of a wavefront are tangent; if the curve degenerates to a point it is a focus|33). This causes only 
a local breakdown of the method and can be handled by a number of methods. At a distance much greater than a 
wavelength away from the caustic the solutions arc still a good approximation to the true solution of the initial value 
problem. 

In contrast, to find asymptotic solutions on a bounded domain D with boundary value conditions, one must 
introduce more than one eikonal at each order in the asymptotic expansion. We will illustrate the important points 
here with Dirichlet boundary conditions on the boundary dD. However the basic argument holds for any linear 
homogeneous boundary conditions and, with minor modifications, for the matching conditions relevant for uniform 
dielectric resonators of index of refraction n with boundary shape dD, within an infinite medium of index n ~ I. For 
the discussion of Dirichlet boundary conditions wc will set the index n = 1 for convenience within the domain D. The 
leading order in the asymptotic expansion of the solution takes the form 

JV 

ij{x) = A,„(a;)e^'^-^-(-) (6) 

m— 1 

with N > 2. It is easily checked that there must be more than one term (eikonal) in the solution in order to have a 
non-trivial solution; if there were only one term in the expression for ^pi^) then any solution which vanished on the 
boundary would vanish identically in D due to the form of the transport equation. 

The question we now address is the following. For what boundary shapes dD in two dimensions do there exist 
approximate solutions of the form Eq. © which are valid everywhere in D except in the neighborhood of caustics 
(which are a set of measure zero)? 

First we note that with Dirichlet boundary conditions we have a hermitian eigenvalue problem and so we know that 
solutions will only exist at a discrete set of real wavevcctors k. In the eikonal theory the quantization condition for k 
arises from the requirement of single-valucdness of ip{x) and will be reviewed bricfiy below. Here our primary goal is 
to show that the existence of eikonal solutions to the boundary value problem is intimately tied to the nature of the 
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ray dynamics within the region D. Moreover for the case of fully chaotic ray dynamics this connection shows that 
eikonal solutions do not exist. We will prove this latter statement by showing a contradiction follows from assuming 
the existence of eikonal solutions in the chaotic case. This argument will be a "physicist's proof without excessive 
attention to full mathematical rigor. 

The proposed solution for ip{x) posits the existence of N scalar functions Sm{x) each of which satisfy the eikonal 
equation, {W Sm{x))'^ = 1 and which, while themselves not single-valued on the domain D, allow the construction 
of single-valued functions 'ip{x) and \7tp{x). Moreover, for the asymptotic expansion to be well-defined, the "rapid 
variation" in 'ip{x) must come from the largeness of fc; i.e. to define a meaningful asymptotic expansion in which 
terms are balanced at each order in k the functions S„ cannot vary too rapidly in space. From the eikonal equation 
itself we know that jV^ml = 1, but we must also have that the curvature V^S*™ <C k for the asymptotic solution to 
be accurate. This condition fails within a wavelength of a caustic, as one can check explicitly, e.g. for the case of a 
circular domain D; but for a solvable case like the circle it holds everywhere else in D. 

It is convenient for our current argument to focus on V?/), instead of ip itself. Consider an arbitrary point xq in D 
where Vipixo) 7^ 0; to leading order in k and away from caustics 



The N unit vectors VSmixo) = p,„ define N directions at xq which are the directions of rays passing through a;o in 
the stationary solution. An important point is that due to the condition on the curvature just noted, these directions 
are constant at least within a neighborhood of linear dimension A = 27r/fc around xq. Choose one of the ray directions, 
call it pi and follow the gradient field V^i to the boundary D. For a medium of uniform index (as we have assumed) 
the vector VSi is strictly constant in both direction and magnitude along a ray. Thus one can find the direction of 
VS*! at the boundary and calculate its "angle of incidence", n ■ VSi, where n is the normal to the boundary at the 
point of intersection. The condition ?/) = on the boundary implies that there is a second term with the eikonal ^2 
in the sum. which satisfies Si ~ S2 and Ai = —A2 on the boundary. As a result tangent derivatives of 6*1,2 on the 
boundary are also equal and together with Eq. this implies that for a non-trivial solution n ■ VS'2 = —n • VSi. In 
other words a ray of the eikonal Si must specularly reflect at the boundary into a ray of another eikonal in the sum, 
which we label 82- Hence we know the direction of VS'2 at the boundary and can follow it until the next "reflection" 
from the boundary. Thus each segment of a ray trajectory corresponds to a direction of V Sm for some m in Eq. Q. A 
ray moving linearly in a domain D and specularly reflecting from the boundary describes exactly the same dynamics 
as a point mass moving on a frictionless "billiard" table with boundary walls of shape dD. Such dynamical billiards 
have been studied since Birkhoff in the 1920's as simple dynamical systems which can and typically do exhibit chaotic 
motion. Thus the problem of predicting the properties of the vector fields Sm is identical to the problem of the 
long-time behavior of dynamical billiards. 

One property of any such bounded dynamical system (independent of whether it displays chaos) is that any 
trajectory starting from a point xg will return to a neighborhood of that point an infinite number of times as t 00 
(the Poincare recurrence theorem[53|). Therefore we are guaranteed that the ray we followed from xq in the direction 
Pi will eventually re-enter the neighborhood of size A around xq. By our previous argument, each linear segment of 
the ray trajectory, corresponds to one of the directions S/Sm and thus when the ray re-enters the neighborhood of 
Xq for Vip to be single valued it is necessary that the ray travel in one of the directions VSmixo) = Pm- There can 
be two categories of ray dynamics: 1) Although the ray enters the neighborhood of xq an infinite number of times it 
only does so in a finite number, N of ray directions. 2) The number of ray directions grows monotonically with time 
and tends to infinity as t —^ 00. We will now show that the general applicability of the eikonal method depends on 
which category of ray motion occurs. 

Let us first consider a billiard dD with fully chaotic dynamics. In the current context "fully chaotic" means that 
for arbitrary choice of xq and the direction pi (except for sets of measure zero, such as unstable periodic orbits) the 
distribution of return directions (momenta) is continuous and isotropic as i — > cxd. Therefore the number of terms 
in an eikonal solution of the form Eq. © would have to be infinite, contradicting our initial assumption that N 
was finite. Thus there do not exist eikonal solutions with finite N for wave equations on domains with fully chaotic 
ray dynamics. A very closely related point was made by Einstein as early as 191 zfl^ (he phrased it as the non- 
existence of a multi- valued vector field defined by the N "sheets" of the functions Sm)- One may ask whether an 
eikonal solution with an infinite number of terms could be defined; this appears unlikely as the amplitudes for the 
wavefronts are bounded below by (A/L)^/^, where L is the typical linear dimension of D, so that only a very special 
phase relationship between terms would allow such a sum to converge. The essential physics of this breakdown of the 
eikonal method is that in a chaotic system wave solutions exist but do not have wavefronts which are straight on a 
scale much larger than a wavelength, hence it is impossible to develop a sensible asymptotic expansion with smooth 
functions 5m- 



N 




(7) 



5 



Returning now to the case of a boundary dD for which the distribution of return momenta is always discrete, this 
means that there exist exactly N ray directions for each point Xq and any choice oipi. In this case the entire spectrum 
of the wave equation on dD can be obtained by an eikonal approximation with N terms of the form Eq. © . The 
quantized values of k are determined by the conditions that the eikonal only advance in phase by an integer multiple 
2tt upon each return to xq and hence the solution is single- valued. The correct quantization condition must take into 
account phase shifts which occurs for rays as they pass caustics. The details of implementing this condition have 
become know as Einstein-Brillouin-Kellcr quantization^^!- 

From modern studies of billiard dynamics we know that both of the cases we have just considered are exceptional. 
The billiards for which eikonal solutions for the entire spectrum is possible are called integrable, and their ray 
dynamics has one global constant of motion for each degree of freedom. For example in the circular billiard both 
angular momentum and energy are conserved and for each choice of Xq and direction pi there are exactly two 
return directions [s^ (see Fig. While the circle is a good and relevant example here, there are other shapes, 
such as rectangles and equilateral triangles for which the method also works; obviously these are shapes of very high 
symmetry. It is also known that an elliptical billiard of any eccentricity is integrable; however this is believed to be 
the only integrable smooth deformation of a cir cle[i|H|. Thus there is a relatively small class of boundaries for which 
eikonal methods work globally; this point does not seem to be widely appreciated in the optics community. 




(a) (b) 

FIG. 1: (a) A typical quasi-periodic ray motion in a circular billiard. The two possible ray return directions for a specific point 
Xo and initial direction pi are shown in red. (b) The Bunimovich stadium, consisting of two semi-circles connected by straight 
segments, for which ray motion is completely chaotic. As the schematic indicates, for any point xo the ray return directions 
are infinite, continuously distributed and isotropic, making an eikonal solution impossible. 

As already noted, the type of boundary shape which generates continuous return distributions for each choice of a;o 
and direction pi correspond to completely chaotic billiards and such shapes arc also quite rare. No smooth boundary 
(i.e. dD for which all derivatives exist) is known to be of this type. A well-known and relevant example for us of such 
a shape is the stadium billiard, consisting of two semi-circular "endcaps" connected by straight sides. Note that the 
generation of continuous return distributions would fail for a point xq between the two straight walls if we chose pi 
perpendicular to the walls generating a (marginally stable) two-bounce periodic orbit passing through xq. However 
this choice represents a set of measure zero of the initial conditions in the phase space. It follows from our above 
arguments that eikonal methods would fail for the entire spectrum in such a billiard (except a set of measure zero in 
the short wavelength limit). 

The generic dynamics of billiards arises when the boundary is smooth but there is not a second global constant of 
motion; this is exemplified by the quadrupolc billiard we study extensively below (see definition in Eq. ijSJ)). Such a 
billiard has "mixed" dynamics; we shall explain what this means and how it is studied in more detail below. For such 
a billiard, depending on the choice of the initial phase space point (xo,Pi), one may get either a finite number N of 
return directions or an infinite number as i — *■ oo. It is not obvious just from our above arguments that this means 
that eikonal methods will fail in such a case. We will skip over this point and simply state that in the case of mixed 
dynamics in principle only a finite fraction of the spectrum could be calculated by eikonal methods. If one can obtain 
a relatively tractable expression for the locally-conserved quantity which leads to a finite number of return directions 
iV, as in the adiabatic approximation of Berry and Robnik[60|. then some progress can be made|49l|. However in 
practice the vector fields VS*™ required are usually too complicated to make such an approach tractable. Thus in 
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practice cikonal methods are not very useful to find solutions of the Helmholtz equation for generic shapes. A related 
but different analytic method, that of Gaussian optics, can be used to calculate a fraction of the spectrum based on 
motion near stable periodic orbits. This method is worked out for dielectric billiards with mixed dynamics in detail 
in Ref. |65j . However both this and the eikonal method fail for a fraction of the spectrum which approaches unity as 
the chaotic fraction of phase space approaches unity. 

Since the traditional analytic methods of optics fail for these systems, what other short wavelength approaches 
exist? The development of short wavelength approximations for mixed and chaotic systems is precisely the problem 
of quantum chaos which has been widely studied in atomic, nuclear, solid-state and mathematical physics over the 
past two decades|3, ll^, IH, H^. Powerful analytic methods have been developed, but with an essentially different 
character than eikonal or Gaussian methods (these techniques are typically referred to as semiclassical methods in 
the quantum chaos literature). The analytic methods in quantum chaos theory are all of a statistical character and 
do not allow one to calculate individual modes. Instead the methods focus on the fluctuating part of the density of 
modes and the statistical properties of the spectrum (e.g. level-spacing distributions). The results are useful in many 
contexts, but less useful in the context of optical resonators and micro-lasers for which a single or small set of modes 
will be selected and one is interested in their emission patterns and Q-values. Therefore it is particularly important 
to develop efhcient numerical methods for calculating the spectrum and modes of such dielectric resonators, and we 
will discuss our method for doing this in section (|VI|) below. 

We are primarily interested in ARC resonators with mixed billiard dynamics as these shapes lead to resonances with 
high Q and directional emission. For such resonators, the ray phase space is not fully chaotic, but is highly structured. 
Moreover the possibility of ray escape decreases the randomizing effect of chaotic motion at t ^ oo. Therefore using 
methods which maintain a connection between the wave solutions and the ray phase space is very helpful. We shall 
describe such a method, known as Husimi projection, in section (jVIIip below. 

III. RAY DYNAMICS FOR GENERIC DIELECTRIC RESONATORS 

Before introducing our numerical method and the Husimi projection method, we review the properties of mixed 
phase space via the surface of section method in the context of the quadrupole billiard/ ARC. This billiard is described 
by the boundary shape: 

R{(j}) = R„{l + ecos2(j)) (8) 

which in the zero deformation limit e = reduces to a circular billiard, which as we have already noted, is integrable. 
Therefore the variation of the parameter e starting from zero induces a transition to chaos. As the perturbation is 
smooth, various results in dynamic a.1 sy stems (collectively known as Kolmogorov-Arnold-Moser theory) imply that 
the transition to chaos is gradual [a. l39l|. The quadrupole billiard displays the typical behavior characteristic of this 
transition. In our initial discussion here we treat the ideal perfectly-reflecting billiard; later wc will discuss the role 
of ray escape in ARCs. 

When the shape is gradually deformed, it quickly becomes unfeasible to capture the types of ensuing ray motion by 
standard ray tracing methods in real space. A standard tool of non-linear dynamics, which proves to be very useful 
in disentangling the dynamical information, is the Poincare surface of section fSOS)|43ll57j. In this two-dimensional 
phase-space representation, the internal ray motion is conveniently parametrized by recording the pair of numbers 
((/)i,sinxi) at each reflection i, where (jji is the polar angle denoting the position of the ith reflection on the boundary 
and sinxi is the corresponding angle of incidence of the ray at that position (see Fig. Ej). Each initial point is 
then evolved in time through the iteration of the SOS map i ^ i + 1, resulting in basically two general classes of 
distributions. If the iteration results in a one-dimensional distribution (an invariant curve), the motion represented 
is regular. On the other hand exploration of a two-dimensional region is the signature of chaotic motion. 

The transition to ray chaos in the quadrupole billiard is illustrated in Fig.|21 At zero deformation the conservation 
of sinx results in straight line trajectories throughout the SOS and we have globally regular motion. These are the 
well-known whispering gallery (WG) orbits for sinx > 1/n. As the deformation is increased (sec Fig. O chaotic 
motion appears (the areas of scattered points in Fig. ^ and a given initial condition explores a larger range of values 
of sinx- Simultaneously, islands of stable motion emerge (closed curves in Fig. O, but there also exist extended 
"KAM curves" [39j (open curves in Fig.O, which describe a deformed WG-like motion close to the perimeter of the 
boundary. These islands and KAM curves cannot be crossed by chaotic trajectories in the SOS. As the transition to 
chaos ensues, a crucial role is played by the periodic orbits (POs), which appear as fixed points of the SOS map. The 
local structure of the islands and chaotic layers can be understood through the periodic orbits which they contain. 
Thus, the center of each island contains a stable fixed point, and close to each stable fixed point the invariant curves 
form a family of rotated ellipses. The Birkhoff fixed point theorem^^ guarantees that each stable fixed point has an 
unstable partner, which resides on the intersection of separatrix curves surrounding the elliptic manifolds. Chaotic 
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Surface of Section for the Qucdrupole at e= 0.12 




FIG. 2: The construction of the surface of section plot. Each reflection from the boundary is represented by a point in the SOS 
recording the angular position of the bounce on the boundary (0) and the angle of incidence with respect to the local outward 
pointing normal (sinx)- For a standary dynamical billiard there is perfect specular reflection and no escape. For "dielectric 
billiards" if sinx > sinxc > 1/n, total internal reflection takes place, but both refraction and reflection according to Fresnel's 
law results when a bounce point (bounce ij^A in the figure) falls below the "critical line" (shown in red) sin^ > sinxc- Note 
that sinx < correspond to clockwise sense of circulation. We do not plot the sinx < region as the SOS has reflection 
symmetry. Below we will plot the SOS for ideal billiards without escape unless we specify otherwise. 

motion sets in at separatrix regions first, and with increasing deformation pervades larger and larger regions of the 
SOS. Akeady at e = 0.1, much of the phase space is chaotic and a typical initial condition in the chaotic sea explores 
a large range of sinx, eventually reversing its sense of rotation. 

IV. FORMULATION OF THE RESONANCE PROBLEM 

A dielectric resonator is significantly different from the closed (Dirichlet) problem due to its openness. In contrast to 
ideal metallic cavities which possess normal modes at discrete real frequencies, dielectric resonators are characterized 
by a discrete set of quasi-bound modesfl^. l4l|. or resonances. As a result, the quasi-bound modes of a resonator are 
characterized by a frequency u = ck and a lifetime r, where c is the speed of light and k = 27r/A is the wavevector 
in vacuum. Experiments on resonators fall into two broad categories, and the presence of quasi-bound modes are 
manifested differently in these two situations. 

In scattering experiments, an incoming field produced by a source in the farfield (spatial infinity) gives rise to an 
outgoing field which represents the response of the resonator, as measured by an ideal detector in the farfield. In the 
ideal case, where absorption is absent, this corresponds to a situation where energy is conserved and hence in this 
situation the EM field has a real frequency, to, which is arbitrary and set by the source. In emission experiments, on the 
other hand, there is no incoming field, but only an outgoing field. As a result, energy is depleted from the system, and 
this process is characterized by decay. The simplest mathematical description of these two experiments correspond 
to the solution of the wave-equation (which is derived from the Maxwell's equations as described in section (|V|l ) 

where the solutions have the separable, time-harmonic dependence 

■f{x,t)^tl;{x)e"^' (10) 

so that ipi^) obeys the Helmholtz equation 



{V^- + n^{x)k^)%l>{x) = {) 



(11) 
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e = 0.05 




FIG. 3: The SOS of a quadrupole at fractional deformations e — 0,0.05,0.11,0.18. The closed curves and the curves crossing 
the SOS represent two types of regular motion, motion near a stable periodic orbit and quasi-periodic motion respectively. 
The regions of scattered points represent chaotic portions of phase space. A single trajectory in this "chaotic component" will 
explore the entire chaotic region. With increasing deformation the chaotic component of the SOS (scattered points) grows with 
respect to regular components and is already dominant at 11 % deformation. Note in (b) the separatrix region associated with 
the two-bounce unstable orbit along the major axis where the transition to chaotic motion sets in first. 



Here, n{x) represents the index of refraction. In general, one can define a complete set of incoming {tpfj, (fc; x)} and 

outgoing modes {ipj^\k;x)} at a given k, in the absence of the resonator. The exact form of these sets is dictated 
by convenience, and in the present discussion we will employ the cylindrical harmonics. 

The two experimental situations at this point are distinguished by two different boundary conditions in the farfield. 
The scattering experiment corresponds to the boundary condition 

i;{x)^ijl-\k;x) + Y,S^,AkWHk;x), \x\ ^ ^ (12) 

and experimentally, it is the scattering matrix 5^1^ (fc), which contains the information measured by the farfield 
detector. In the typical case Sfj,v{k) will display sharp peaks at a discrete set of real wavevectors ki (in case of isolated 
resonances; see the back panel (Re [kR]-I plane) on Fig.^. This is the signature of long-lived quasi-bound modes with 
frequency uj = ckf, their lifetimes r are encoded in the functional form of the peaks, which in general is of Fano shape, 
with direction-dependent parameters. This makes scattering boundary conditions less convenient for the extraction 
of the quasi-bound mode structure. 

The emission experiments are modeled by the outgoing wave boundary conditions at infinity 

V'«(a;)~E>(^0^^^^(^-=^)' \x\^^ (13) 
1/ 

This form at infinity docs not permit solution for any real k as it manifestly violates current conservation. Instead 
the solutions of Eq. I|ll|) . indexed by i, exist only at discrete complex wavevectors ki — Ki -\- iVi. The connection to 
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quasi-bound modes is then direct; the real part gives the quasi-bound mode frequency u>i = cKi and the imaginary 
part represents the Ufetime of the mode, = l/cFj. Here, we will use the radiation boundary conditions exclusively, 
and quote the dimensionless complex variable kR instead, where R is the mean radius of the resonator. In an active 
medium, one may think of these resonances as being pulled up to real wave-vectors by the gain. Although the 
description of actual (stationary) laser modes requires the solution of a non-linear wave-eauation[6^ . we will focus 
here on the problem of linear resonances, an approximation which is often used in laser theory 




Re[jtR] 



lm[jtR] 



FIG. 4: A comparison of scattering and emission pictures for quasi-bound modes. Variation of the intensity scattered off 
a dielectric circular cylinder with the wavenumber k of an incoming plane-wave is plotted on the back panel (7 — Re [kR] 
plane). The intensity is observed at 170° with respect to the incoming wave direction; it would look significantly different in 
another direction due to interference with the incident beam. The complex quasi-bound mode frequencies are plotted on the 
Re [kR] — Im [kR] plane. Notice that the most prominent peaks in scattering intensity are found at the values of k where a 
quasi-bound mode frequency is closest to the real-axis. These are the long-lived resonances of the cavity. Also visible is the 
contribution of resonances with shorter lifetimes (higher values of Im [kR]) to broader peaks and the scattering background. 

The relation between the linear emission and scattering picture is easily visualized in the extended complex wavevec- 
tor space of the scattering matrix S^v{k), depicted in Fig. 0] The discrete quasi-bound wavevectors hi are the poles 
of Sfi^{k). As can be seen from the figure, in general there are multiple quasi-bound modes contributing to a given 
resonance peak, but the quasi-bound modes which are closest to the real-axis lead to the sharpest peaks (some of 
which might not even be resolved in the scattering profile). Note that via Eq. (|10|l . the quasi-bound mode solutions 
damp in time. An important experimental value often quoted is the Q-valuc of a resonator, which is defined by the 
number of cycles of the optical field at frequency lo to decay to half of its value, and thus can be related to quasi-bound 
mode parameters by the relation Q = wr = — 2Re [kR\ /|Im [kR] \ . 

It is possible to generalize eikonal theory to calculate the quasi-bound states of complex k of dielectric cavities 
of integrable shape[6^. The dielectric boundary conditions then reduce to ray trajectories which still propagate on 
straight lines and undergo specular reflection, much like in the case of billiards. The additional feature is that dielectric 
billiards exhibit ray splitting at the boundary, and give rise to both a refracted and reflected ray with amplitude and 
direction obtained from the application of the laws of Snell and Fresnel for a fiat dielectric boundary. The transport 
equations for the amplitude have to be supplemented by an additional complex multiplicative factor at each encounter 
with the boundary. The practical implication is that the ray motion as displayed on the SOS can be used for the 
dielectric problem, when augmented with an escape condition. The escape probability will be exponentially small 
in wavenumber k for angles of incidence above the critical angle Xc = sin~^ 1/"-, since its due to a tunncling-like 
process [43. l64|. This condition is demarcated by the line sinx = sin^c h^ the SOS; any ray falling below this line 
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refracts out with the probability given by the local Fresnel law of refraction (assuming a TM mode) : 



T(sinx)^ .^^7^"? (14) 

Vl-sm X + V sm Xc-sm x 

providing a classical loss mechanism and leading to finite lifetimes of the corresponding modes. Thus the openness 
of dielectric resonators is not the cause of the failure of standard methods; these methods fail for chaotic resonator 
shapes for the same reasons of dynamical complexity discussed in section ^1 for the Dirichlet case. We now discuss 
the reduction of the wave equation for dielectric resonators to the Helmholtz equation to lay the groundwork for the 
numerical method of section IVTI which allows us to solve for the resonances of chaotic shapes. 



V. REDUCTION OF MAXWELL'S EQUATIONS 



Consider the problem of excitation of electromagnetic waves in an infinite dielectric rod of arbitrary cross-section 
(see schematics in Fig.lSJ, which is extended along the z-axis. In practical situations, the structure is of finite extent 
and there are planar end-caps which makes it truly a resonator. In other cases, it's a fiber-optic cable of practically 
infinite extent. In any case, we will for now assume translational symmetry along z-axis, and we will show later that 
this is a perfectly valid assumption for the modes of relevance to us. 




FIG. 5: Illustration of the reduction of the Maxwell equation for an infinite dielectric rod of general cross-section to the 2D 
Helmholtz equation for the TM case (E field parallel to axis) and fcy = — 0. 



Assuming harmonic time-dependence of the fields and no surface currents and charges. Maxwell's equations for the 
system yield the Helmholtz equations: 

{V' + n^{x)k^)(^^^^0 (15) 

where A: = 27r/A = w/c is the wavevector in vacuum; n ~ y/JIe is the index of refraction, ^ is the permeability and 
e is the dielectric constant of the medium, which are in general function of position. We will assume /i — 1, so that 
= e. 

The translational symmetry along the z-axis allows us to express the .j-variation of the fields as 

E{x) ^ E{x,y)c'"'''' (16) 

Following Jackson|32j. we separate the fields and operators into components parallel and transverse to the z-axis and 
write out the transverse projection of curl equations: 

ink^E±+ikz X B± = V±E^ (17) 
ink^B±_ - inkz X E± = Vx^z (18) 

It's evident from these four (scalar) equations that Ez and Bz are the fundamental fields we should be after, and 
that once they are determined we can solve for Ej_ and B±. Thus, the Maxwell's equations themselves completely 
decouple, which was already obvious from Eq. H15() . The actual complication of solving the vector Helmholtz equation 
stems from the fact that the boundary conditions are coupled. The Maxwell boundary conditions are 

Ox {El- E2) =0, u- {njEi - nlE2) = (19) 
i> X (Bi - B2) = 0, 7>-(Bi-B2)=0 (20) 
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in the absence of surface currents and charges and for a hnear, isotropic medium. The subscripts denote the media 
on respective sides of the interface. i> is the unit normal on the interface, pointing towards out from the cyhnder. We 
will assume n\ ~ n > n-i ~ \. Note that these are six conditions altogether. Focusing on the scalar fields E^, Bz, we 
have from the equations involving the cross-product with i), Eiz = E2z and Biz = B2z- Another pair of boundary 
conditions can be found by projecting Eq. ((T5|) onto the diad {0, s) defined on the boundary of the cross-section dD 



dEi^_dE2z_ ^ kz f 1 dBiz dB2z 
ds ds k \n dv dv 

dEiz dE2z kz f 1 dBiz dB2z 



dv dv k \n ds ds 



(21) 
(22) 



We are interested in the long-lived modes of the resonator. Modes with a finite kz correspond in short wavelength 
limit to rays which spiral up and down along the cylinder walls and escape through the end-caps by refracting 
out. Thus under most circumstances the longest lived modes have kz ~ 0, and correspond to modes which are 
effectively two-dimensional, i.e. can be expressed by dynamics of rays on the cross-sectional plane. In that case, the 
boundary conditions Eq. also become diagonal and we have a complete decoupling. We will choose to work with 
tjji{x,y) = Eiz{x,y), corresponding to TM polarized fields, for which the problem reduces to the two-dimensional 
Helmholtz equation for the scalar field ip with continuity conditions 

(Vi+n2fc2)V^,;(a;,y) =0 (23) 

Vi\dD = ^'2 an, an = -^Iod (24) 

dv dv 

Note that this boundary value problem is equivalent to that of the stationary Schrodinger equation of quantum 
mechanics. Hereafter, we will drop all references to the original three-dimensional and vector character of the problem 
and work with Eqs. 

VI. SCATTERING QUANTIZATION-PHILOSOPHY AND METHODOLOGY 

In this section we will describe a numerical method to solve Eq. H24|) , which is both efficient and physically appealing. 
Our approach is a generalization to open systems (specifically, dielectric resonators) of the scattering quantization 
approach to quantum billiards |l5lll6l| . This approach is based on the observation that every quantum billiard interior 
problem (Helmholtz equation for a bounded region with Dirichlet/Neumann boundary conditions) can be viewed as 
a scattering problem, and the spectrum can be uniquely deduced from the knowledge of the corresponding scattering 
operator. In the case of closed systems, the internal scattering problem can be mapped rigorously to an external 
scattering problem |l7l|. and the resulting (exact) scattering matrix is unitary. For the dielectric resonator problem 
with radiation boundary conditions, we will see that the corresponding scattering operator is inherently non-unitary, 
reflecting the physical fact that we are dealing with a leaky system. Thus we will define below a new "S-matrix" 
which is non-unitary and distinct from the true S-matrix describing external scattering from the system. We retain the 
terminology "S-matrix" nonetheless because of the conceptual similarity to the quantum billiard method of [Tsl Il6l | . 
The generalization of this approach to dielectric billiards was flrst made in Ref.^3|, however without the efficient 
algorithm presented below. 

We assume that the resonator is bounded by the interface dD of the form r = R{4'), where R{4>) is some smooth 
deformation of the boundary such that there exists only one point of the boundary for each angle cf). We decompose 
the internal and external fields into cylindrical harmonics with a constant k 

oo 

V'i(r,0) = (amH+(nfcr)-t-/3™H-(nfcr))e™'^ ^ < (25) 

m— — 00 
00 

^2(r,0) = (7™H+(fcr)-t-<5,„H-(A;r))e™* ^ > ^(^) (26) 



m— — CXD 



Each of the terms 

V'±(r,(/))=H±(nA;r)e""^ (27) 
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in the sum is a solution of the appropriate (interior or exterior) Helmholtz equation, but does not satisfy the matching 
conditions by itself. Note that {V'm} forms a normal basis in the infinite space. Owing to the completeness of this 
basis, the expansion is exact for r < Rmin and r > Rmax as long as the sum runs over an infinite number of terms, 
where Rmin and Rmax are the lower and upper bounds of R{4>) respectively. The assumption that the expansions 
can be analytically continued to the region Rmin < r < Rmax is known as the Rayleigh hypothesis [s^. It has been 
shown that for a family of deformations parametrized by e, there is typically a critical deformation Cc, beyond 
which the hypothesis breaks down because the expansion ceases to be analytic in the region Rmin < r < Rmax- For 
the deformations Eq. (jS)), this happens long after the shape becomes concave; wc are not interested in this regime. 
Although this issue seems thus to be resolved, we shall see that precursors of the non-convergence emerge in the form 
of numerical instabilities for e < e^- 

We will assume that Sm ~ (no incoming waves), thus confining our attention to quasi-bound modes. Turning 
to the interior expansion, the regularity of the solution at the origin requires that we take am = Pm, but we will 
not implement this condition at this stage. The continuity conditions Eq. (|24() give us further relations among the 
remaining coefficients: 

= M<f>.Rm (28) 
i0,_R(0) - \4>M4>) i^y) 

In Eq. H29|) . we have replaced the normal derivative condition by the radial derivative condition, because Eq. H28|) 
shows that the tangential derivatives are also continuous. Note that this latter set of equations containing radial 
derivatives is equivalent to the set of equations (|24|l using normal derivatives. 
These conditions can be written out as 

OO CXD 

^ (a,„H+(nfci?(0))-HA„H-(nfci?(</))))e™^ = ^ 7mH+ (fci?(0))e™^ (30) 

m— — OO rn— — OO 

OO OO 

^ (a™H+'(nfci?(0))-f/3,„H-'(nfci?(</>)))e™* = J2 7™H+'(fci?(</>))e™^ (31) 



n 

m— — OO 



We multiply both sides by Wn{<t>)c '""^ and integrate with respect to (j) to get a matrix equation for the coefficient 
vectors \a), and I7) 

n+\a)+HY\l3) = W+I7) (32) 
Vnt\a) +VH^\0) = -Vn+h) (33) 



Various choices of the weight function w{(f)) arc possible^3|; here we choose ?«(0) = 1. The matrices in Eq. H32I I33|l 
are defined by 

[Hf]^^ = r dcf,lltjn,kR{cf,)y^"^~''>^ (34) 





r2TT 



[Vnf]^^^ - I d4>lit:{n,kRi4>)y<^"^-'^^ (35) 

Eliminating I7) between Eq. and Eq. we obtain 

S{k)\a) = 1/3) (36) 

where the matrix S{k) is given by 

s{k) ^ [nivniy^vni - [nty^n^] [{nt)~^nt - n{vnt)~^vnt] (37) 

As noted earlier, this S-matrix is different from the standard external scattering matrix introduced in Eq. I|12l) . It is 
straightforward to check that, for real k, S{k) is non-unitary. Consider now the eigenvalue problem of S{k) 

Sik)\a) =c"^\a), (38) 

where for real k the phase (p is complex. Once we find a complex kg where one (or several) of the (p is a multiple of 2tt, 
we have \a) = which is exactly the condition of regularity at the origin. This is the quantization condition which 
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FIG. 6: Schematics describing the quantum Poincare mapping induced by the internal scattering operator, see discussion in 
text. 



will provide us with the quantized eigenvalues and eigenvectors (kq, \a^'^')) which allow us to construct the resonant 
solutions of the interior and exterior problem we set out to find. This condition is often expressed in terms of the 
secular function C(fc)01 given by 



C(/fc) =dct[l-5(fc)] 



(39) 



The spectrum is obtained as the zeros of the secular equation ^(fc) = 0. As noted, the values kq for which we obtain 
a unit eigenvalue and the secular function Eq. (|39|l has a root, is always complex and the eigenvalues of S{k) arc not 
pure phases, ip € C. The practical upshot of this is that this requires a two-dimensional root-search for the equation 
(^{k) — 0. An often employed numerical procedure involves a sweep in the complex /c-planc of the singular values of the 
operator T{k) = 1 — S{k), with proper care of the numerical null-space of T(A:)0. This requires several calculations 
of the entries of T(fc) and its singular value decomposition per quantized state. In the next section, we will represent 
an efficient root-finding method, which ideally requires two diagonalizations per nkR quantized states. Before doing 
that however, it's worthwhile to investigate the structure of S{k) based on simple physical considerations. 

A physical interpretation of the internal scattering operator S{k) and its eigenvectors can be given even off- 
quantization {f{k) ^ 27r') [T^l3^ . We can visualize this approach in our case by dividing the interior of the resonator 

into two subdomains joined along the curve C, which we take to be circle of radius Rc ^ Rmin, and considering 
it as a boundary at the junction of two back-to-back scattering systems. We furthermore introduce a tiny metallic 
inclusion of radius 6 at the origin (this is introduced for the sake of the argument and can be omitted) . This is our 
first scattering system, which scatters an incoming wave |/?) into |a) via the scattering operator Ss 



Ss{k)\a) 



(40) 



Ss{k) is exactly the exterior scattering operator for a metallic circle (immersed in a medium with index of refraction 
n): 



[Ss{k)]r 



n^jnkS) ^ 



(41) 



The second scattering system is the boundary itself, scattering an incoming wave (with respect to the boundary) 
\a) into 1/3), and the scattering operator for this system is simply S{k) whose form is given in Eq. l|S7|) . Consider 
now a whole cycle, starting with the state \a) on C, being first scattered off the tiny circle, then from the boundary 
returning to C again (see Fig. E)). The resulting scattered vector is S ■ Ss\a). Now, as k6 — > 0, we have Ss 1, 
and the resulting scattered vector is S\a). Because the individual normal modes V'm ™ our expansion correspond to 
ray trajectories which have a well-defined angular momentum sinx = nkkc ' ^^^^ mapping S\a) can be interpreted 
as a wave analogue of the Poincare SOS mapping on the section C, parametrized by ((/>, sin^). This link has been 
fruitfully used to obtain short wavelength forms of the scattering operator iS(fc), for various closed systems js^- We 
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FIG. 7: A gray-scale representation of the scattering matrix Eq. 13711 . calculated for a quadrupolar resonator at e = 0.1 
deformation, n = 2.5 and nkR = 40. The number of evanescent channels used in the calculation is Aev = 15. Note the strong 
diagonal form for |m| > nkR. The spread around the diagonal is proportional to the deformation. Here the internal scattering 
couples approximately 20 angular momentum modes. 

will not pursue this approach here, but will make use of this visualization to develop a meaningful truncation scheme 
for a numerical implementation of our method. 

First of all, at a given k, an angular momentum eigenstate ipm: for which m > nkRmax is a closed channel for 
the internal scattering system, because it corresponds to classical motion with a circular caustic of radius larger 
than Rmax- Such channels are called evanescent^ and are not not expected to be scattered significantly. In fact, 
a plot of the matrix S{k) in Fig. [3 reveals that as m grows beyond a critical value rUc ~ nkRmax, the scattering 
matrix becomes strongly diagonal i.e. [5(fc)]„„„' « 5„„n' for \'m\^ |m'| > me- Furthermore, there is a transition region 
nkRmin < I'T^I, |m'| < nkRmax, where the matrix is heading towards diagonality, and this region corresponds to 
evanescent components which undergo an enhanced scattering because they overlap significantly with only certain 
regions of the resonator. This region grows with the deformation of the resonator, and consequently, Aei, evanecent 
channels have to be included in the number A of (positive) channels contributing to a given internal scattering matrix. 
Deonoting the critical matrix size at the evanescent channel boundary by Kgc — \nkRmin\ (!■] stands for the integer 
part), the size of the 5'-matrix is then Ntmnc = 2A + 1, with A = A^c + ^ev 

VII. ROOT-SEARCH STRATEGY 

A typical run at nkRo = 106 for e ~ 0.1 produces the eigenvalue distribution {e"^''} plotted in Fig. |S1 in the 
complex z = e*"^ plane. We will denote tp = 9 + ir], where 9 and rj are real numbers, so that \z\ = exp(— r;). Note that 
= InkRo^l — e)] = 93 and we have included Agy = 55 evanescent channels. The handling of numerical stability 
issues relating to the inclusion of such a large number of evanescent channels is outlined in section 

Our first observation is that all the eigenvalues are strictly distributed within the unit circle \z\ = 1, i.e. Im [ip] < 0. 
This is because of the restriction of solutions to outgoing waves only. Furthermore, there is an accumulation of 
eigenvalues on the boundary of the circle, particularly at 9 = 2tt^ . As we have established, an eigenvalue for which 
(k) = 2tx within a given numerical precision yields a quantized mode of the resonator. However, we should resist 
the temptation to simply take all the scattering eigenstates whose eigenphases are (/s w 27r to be quantized. As 
was pointed out in Ref.[la| in the case of a closed system, there is an accumulation of scattering eigenphases at 
Lp K, 27r"'", which do not correspond to proper physical eigenmodes of the resonator. These are modes, which are 
primarily composed of evanescent channels, and can easily be distinguished from regular modes, because of their lack 
of k-dependence, as we shall see below. 



15 



oo O 

^ o o o 

o ^ o 

— 


n 

°o °°o \ 

° <> 



OoO . 

'6 

w ° ° % 

°°o ° ® 1 


So® o 
• Co 

I ° °s 

\ ° o 

•,0 off n 

„° c 

\ 

01.,. 


° °" s i 
i 

OS 



0*0 -° / 

n -* 

1 1 







FIG. 8: Distribution of scattering eigenvalues (red circles) in the complex plane for nkR — 106, e — 0.12, n = 2.65. Blue dashed 
line is the unit circle \z\ = 1. Long-lived states have the modulus of the eigenvalue very close to unity, i.e. the eigenphase has 
only a small imaginary part rj. 



A. Zero deformation-Case of the rotationally symmetric dielectric 



A lot can be learned by way of a simple example. We will consider a case where we know the exact solutions, 
namely the dielectric circle. The exact eigenstates of the scattering matrix for the circle can be given a precise 
physical meaning in terms of classical processes in the short wavelength limit. They correspond to motion with a 
conserved angular momentum, or in terms of our notation in section , a given impact angle sin x on the dielectric 
interface. The resulting scattering matrix is diagonal in the angular momentum representation. This signifies the 
fact that a "channel" with a given m upon encountering the boundary will be scattered to the same channel to, 
corresponding to specular reflection. The scattering matrix can be written as 

[5(fc)]mm' = ~Sjnm' ^''^^^^J^, frnik) (42) 



where the function fm{k) is given by the following expression 



fm{k) 



1 — n 



H+'(nfci?) H+ (fci?) 



H+ (nfci?) H+'(fci?) 



1 - n 



H7;(nfcfl) H+(fci?) 
H- (nfci?) H+'(fci?) 



(43) 



This form in terms of the particular ratios of Hankel functions will help us simplify the expressions considerably in 
the asymptotic limit nkR — > oo. Notice that when fm{k) ~ 1, 

[Sc{k)Un' = -^f^ri^s,^^, (44) 

H™(nfci?) 

is the external scattering matrix for the closed circular cavity, which is unitary. Then our quantization condition 
[Scik)]m.m,' = 1 yields 

JrninkR) = (45) 

which is the exact quantization condition for wavevectors nk of a metallic cavity. Hence in the form Eq. (|42(l . the 
corrections due to the openness of the system are lumped into the factor fm{k). 

Let's first consider the diagonal elements of Eq. I|42(l for to > nkR. We will use the notation a = cosh" ^(m/nfci?), 
a' = coshr ^{m/kR), (3 = cos~^{m/nkR) and 5' = cos~^(rn/kR). Note that a' > a ^ 1. Using the large-order 
asymptotic representations for Bessel functions |]J and with proper attention on exponentially small terms, it can be 



shown thatl 

[S{k)Um ~ 1 + i(l + 2n)e-2-" (46) 
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for m ^ nkR. As noted, these entries correspond to scattering of evanescent channels and result in eigenphases 
exponentially close to zero, (p ^ {1 + 2n)e~^™". Thus, the accumulation of eigenphases on the unit circle close to the 
quantization point = 27r in Fig.|Hlcan be linked to such extremely evanescent channels, which arc not the physical 
modes of the cavity. These modes can be interpreted as creeping waves, which are evanescent modes which cling to 
the surface of the resonatorjs^- Note that the number of such scattering eigenstates depends strongly on our choice 
of Aev in our numerical implementation. 

Next, wc will look at the internally reflected channels. These are obtained for the entries kR < m < nkR. The 
asymptotic form of the corresponding matrix elements are|64j 



1 - i 2?i sin /?e-" 



-2ma 



n sin/5 



(47) 



where Q, which is identical to the closed-circle eigenphase, is real and given by 



e(fc) -2m(/3-tan/3) 



(48) 



These channels yield eigenvalues which accumulate exponentially close to the unit circle \z\ = 1, but unlike the 
evanescent modes Eq. (|46|l . with arbitrary phases. Note that the exponentially small difference from |z| = 1 represents 
the evanescent leakage which vanishes in the short wavelength limit. 
It's possible to assign a velocity to these eigenphases in fc-space: 



de 



d{nkR) 



2sin/3 + 



ikR 



> 



(49) 



A useful observation at this point is that this velocity is twice the cosine of the conserved ray impact angle x = ^/2 — /? 
in the circular billiard corresponding to the motion with angular momentum m (see Fig.EJ. 




FIG. 9: Geometric representation of the angle /3; the velocity of the eigenvalues in the complex plane is 2 sin /3, which is also the 
chord-length of the corresponding ray. We note that for the diametral two-bounce orbit the speed is maximal (corresponding 
to the minimum free spectral range) while for whispering gallery modes the chord length is minimal and the free spectral range 
is the largest. 



The picture this entails is the following: When we slowly increase fc, the individual eigenphases move with an 
approximately constant but mode-dependent speed given by Eq. H49|l counter-clockwise around the unit circle. Each 
time one of the eigenphases passes through ip = 27r, the quantization condition is fulfilled and the resulting eigenvector 
is a quantized mode of the resonator. Hence, the eigenvectors of S{k) can be assigned a physical meaning and identity 
even when k is not tuned to resonance ip{k) = 2tt. In the present case, they correspond to totally internal reflected 
whispering gallery modes. 

Last, we investigate the classically open channels, which corresponds to rays which are refracted out. In this regime 
m < kR and 



sm p' + n sm p 



(50) 



Thus, 



Note that the algebraic prefactor is the Fresnel reflection factor for a ray coming in at an angle Xi = f ~ 
the proximity of the scattering eigenphase to the unit-circle is a measure of the lifetime. The smaller the radius of 



17 



the cigcnphase, the smaller is the associated lifetime. As we change k, the variation of the eigenphase of a given 
solution will be dominated by the phase- factor e'® . The path to quantization goes thus by first increasing Re [k] until 
Re [9] ~ 2tt, and then adding a small imaginary part iAk so that 



|g-«(fc+jAfe)| ^ sin/3'-nsin/? 

sin f3' + n sin (3 



(51) 



driving the eigenphase right to the quantization point. From this condition, we can extract an approximate value for 
the imaginary part of the quasi-bound mode which will result: 



Im \nkR] 



1 



2sin/3 



log 



sin f3' — n sin (3 



sin P' + n sin /? 



(52) 



This is precisely the lifetime of refractive WG modes due to Fresnel scattering, which can be obtained using different 
methods mil. 

The crucial point here is that these statements are only valid for an interval of the order of a mean-level spacing, 
so that (3 is approximately constant 



d(3 



d{nkR) 



O 



1 



(53) 



Furthermore, the assumption that Im [nkR] <§; Re [nkR\ is also implicit in these derivations. These procedures have to 
be implemented carefully because of the Stokes phenomenon |^|^ in the asymptotic expansion of the Hankel functions 
with complex argument. However, as long as the latter condition is satisfied, these estimates are valid. 



B. Deformed dielectric resonators 



In light of our findings for the undeformed case, it is possible to develop a powerful search strategy for the general, 
deformed case. The reason behind our ability to "track" the scattering eigenphases through quantization in the case 
of the circular resonator was the fact that the angular momentum channels didn't mix when we changed fc, owing 
to the diagonality of the scattering matrix over all k i.e. there we had a good label m which was conserved. This 
will not be the case when we deform the resonator. For small deformations, the internal scattering matrix S{k) will 
remain approximately diagonal, with fluctuations due to intcr-channel scattering. The resulting eigenstates will show 
a broadening in their angular momentum distributions. In that case, one can still define an average phase velocity 
given by 

'^^ 2sinP (54) 



d{nkR) 

defined by the average angular momentum in 



A 



^ = ^°^";^ -=^A^EH«™I' (55) 

-A 

At first sight, there is no reason for such a solution to persist over a given interval A/c. Following Ref.^^, we suggest 
that the scattering eigenvectors have an identity beyond a given fc-value, and more importantly, that the resonances, 
the quantized modes, have an identity even when they don't fully satisfy the boundary conditions. We can quantify 
this statement by defining a simple scalar product between eigenvectors of the internal S-matrix at different k: 

(a(fc)|a(fc + Afc)) =^a„(fcX(fc-f Afc) (56) 

m 

Then our claim is tantamount to the adiabaticity of {a{k)\a{k -I- A/;;)). The reason this is possible lies in the subtle 
correlations among the matrix elements induced by the underlying classical motion in the short wavelength limit. We 
have already emphasized the connection between the scattering matrix in the short wavelength limit and the classical 
SOS map. As long as there are invariant curves in the SOS, which we have seen is guaranteed by the KAM scenario 
for near-integrable deformations in section (|III|I , there will be eigenstates of the scattering matrix which will display 
the aforementioned adiabatic behavior. 
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FIG. 10: The overlap calculated for a set of states in the interval nkR = 106 — 107.5, for e = 0.12 and n = 2.65. The associated 
classical structures are found from the Husimi projections of the respective states (see Fig. IIH . The schematics identify the 
ray orbits with which they are associated (see discussion in text). At this deformation the short two-bounce orbit is stable, the 
long one unstable; the diamond orbit is stable and the fish and triangle orbits are unstable. 



In Fig. 1101 we trace the overlap Eq. (|56|l of a set of eigenvectors in an interval of the order of a mean- level spacing. 
First, a diagonalization of S{k) is performed at a fco, the eigenvectors determined, and then further diagonalizations 
are performed at regular intervalls k = ko + jAk, where nAkR = 0.03. At each step, there is in general a single state 
having markedly higher overlap with the respective original state at fcg than the others and that value is plotted. The 
result shows that an adiabatic identity can be in fact defined for certain states. This procedure allows the tracking of 
majority of the states, as long as the deformation is not too large. In fact, it's possible to show that 

{a{k)\a{k+jAk)) = l+jnAkR-0{^) (57) 

At this point it may be helpful to clarify what we mean by the "identity" of a state in the chaotic case in which 
the state is not associated with a stable periodic orbit or a family of quasi-periodic orbits. In Fig. 1111 we show both 
a real-space solution for the electric field of a TM resonance and its projection onto the surface of section using the 
Husimi-SOS projection technique defined in section ljVIII|) . This method allows one to associate any solution, even 
a non-quantized one, with a region in the SOS and hence with an approximate ray-dynamical (classical) meaning. 
Furthermore, at the values of nkR at which we work, often these chaotic states are associated with unstable periodic 
orbits or their unstable manifolds (this is the case in Figs. 1101 and lll|) : in Fig. ^2 the top states are associated with 
the unstable fish orbit and the bottom states are associated with the unstable manifolds of the unstable two-bounce 
orbit along the major axis of the resonator. States localized on unstable periodic orbits have been termed "scars" in 
the quantum chaos literature and will be discussed further in section (|X|| . 

It turns out that one can extend this strategy to higher deformations, where the SOS displays large chaotic com- 
ponents, with proper attention to eigenstates which have an appreciable overlap with chaotic regions. A typical 
scenario which is encountered is the avoided crossing of two scattering eigenvectors. This is captured in Fig. 1111 
where two eigenvectors are traced over a mean-level spacing. Originally, the two states are well-distinguished; they 
have approximately zero overlap with each other. They have different classical meaning as well as shown by their 
Husimi projections on the SOS (see section (|VIII|) for definition). One state is associated with the border of the 
stable bouncing ball region of the SOS and has no intensity near <^ = 0, tt; the other is concentrated in the separatrix 
region associated with the unstable period two orbit along the major axis (we have plotted its unstable manifold for 
reference). At the crossing they perturb each other strongly, and an approximate superposition state results. However, 
if we continue changing k, the states emerging from the avoided crossing will still have a pronounced overlap with the 
states before the crossing. Notice that the overlaps are calculated with reference to one of the original states \ao). 
This example represents a case where a numerical tracing algorithm has to be properly conditioned. 
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FIG. 11: Two eigenvectors are traced by the criterion that the overlap is largest in two consecutive iterations. The figure shows 
the overlap of the two sets of states resulting with respect to one of the initial states, |ao). Away from the avoided crossing 
the states have distinct classical meaning as discussed in the text; they exchange "identity" at the avoided crossing. Both the 
real-space electric field intensities are shown (false color scale) and the Husimi-SOS projections of the states before and after 
the avoided crossing. 



After having established that we can assign an identity to the scattering eigenvectors as k varies, we next investigate 
how precisely the corresponding eigenvalues move within the complex unit circle as we vary k, both through real and 
imaginary values. Fig. I12b 'l shows such a tracing of several representative states. First, the initial eigenvalues are 
followed while varying the real part of fc; each of the eigenvalues follow approximately a circular trajectory, followed 
by a pure imaginary change in k resulting in the eigenvalues following an almost precisely radial path. We write the 
radius of the complex eigenvalue as |e*''''-'^^| = e''^'^-' and call 9 the angle in the complex plane for the eigenvalue e*^'*^^. 

This simple behavior can be understood from the fact that the classical channels (of angular momentum in our 
case) in the expansion preserve their identity over a mean level spacing, and the weight of these channels embodied 
in the expansion coefficients am change only O (^jjr^^- In conclusion, the radial and angular speeds of the eigenvalue 
are approximately "decoupled". This speed is to high accuracy constant for the eigenphases, i.e. for the log of the 
eigenvalues as shown in Fig. 112b . c. 

We have developed an efficient numerical algorithm to determine the quasi-normal modes of an smoothly deformed 
dielectric resonator based on all of these observations: 

1. A diagonalization of S{k) is performed at a given fc, and Ntmnc eigenphases and eigenvectors are determined, 
denoted by |ao''), i = 1, . . . Ntmnc] (»7i|ao ^) = a™ . 

2. A second diagonalization is performed aX k + A/c, where Afc is a small complex number so that |Afc| <C k. 

3. Approximate radial and angular eigenphase speeds are determined. 

4. Assuming the constancy of the individual speeds, an approximate quantization wavevector kq^ is determined 
for each of the initial eigenvectors. 

5. Finally, the quasi-bound modes are constructed by 



A 

i^t\r,cl,)= a^)J™K^V)e™'^ (58) 

r)i=— A 

Due to the small change in {am} with k, we simply use their non-quantized values in the expansion with the 
extrapolated fc- value kq \ We have checked that it's important to use kq '^ instead of the original k. 

In the ideal case this means that Ntmnc nkR quasi-bound modes are found in only two diagonalizations. In 
practice, this ideal limit is not fully attained. But. depending on the deformation and the value of nkR, a large 
fraction of the quasi-bound modes can be calculated approximately in this manner. Table Fig. shows a typical 
run and the quality of the results compared to "exact" solutions. For increased numerical stability we have found it 
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FIG. 12: Several representative eigenvalues z (corresponding to states associated with the ray orbits above and below via color 
code) traced in the complex plane as one changes the real and imaginary values of k. First Re [fc] is varied resulting in the 
circular arcs of fixed radius (Im[2:]); subsequently Im [fc] is varied resulting in a radial motion and fixed Re[fc]. On the right 
hand side we show the constancy of the derivative of the phase angle Q with respect to Re [nfcJ?] and of the derivative of the 
logarithm of the radius r\[k') with respect to Im[nfc_R], implying constant speed of the eigenphases as a function of k in the 
complex plane. The simulations are performed at e = 0.12 quadrupolar deformation, n = 2.65. 

convenient to use the same algorithm in the context of solving the equivalent generalized eigenvalue problem for the 
system; this is discussed briefly in the Appendix IXI 

The implementation of the algorithm can be adapted to the particular result of interest. In fact, when the quan- 
tization of a single state is desired, a more exact eigenphase quantization can be performed by multiple-scans with 
update of the speeds, quite like Newton's root-search method. 

As already hinted, the method is most powerful when applied to calculate classically meaningful quantities (such 
as the Husimi projection, see section (|VIIII| below) for which finding the quantized k values is not important. This 
is most evident in considering the physical observables which are unique to the open systems, the farfield emission 
patterns and the boundary image fields. The calculation of these observables will be described below. 

On the other hand, it has to be pointed out that an attempt to plot an internal solution away from quantization is 
not meaningful because of the nature of the Hankel function basis. The internal field would be 

= ^ a„(H+ (Tifcr) + e^'^H" (nfcr))e^™'^ (59) 

m 

The existence of the factor e!'^ exposes the z"™ singularity at the origin due to the Neumann components N-m(z) = 
i(H~ (z) — H~(z)). That is one reason that we must use kq in Eq. (|58|l . This limitation is an artifact of the particular 
basis used, and does not arise for example for the case of an expansion in cartesian modes. 
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TABLE I: A typical run at kRo = 40, e = 0.12 and n = 2.65. The first column represents the predicted value assuming a 
constant speed as determined from two successive diagonalizations separated by AkR = 10~^ + ilO~*. The second column is the 
error of this prediction, obtained by a full root search, measured by the distance in the complex plane between the eigenphase 
and the quantization point. The last two columns are the overlaps of the original eigenvectors (internal and external) with the 
actual quantized ones. 







kR 


9 






le'^ - 1| 


{aq\ao) 


(79l7o> 


40 


530139923096 


~iO 


128341300297e 


- 


03 





1049356E- 


01 





935212779 





815777069 


40 


354640960693 


-iO 


346842617728e 


- 


02 





2512625E- 


01 





800668216 





888131129 


40 


362663269043 


-iO 


262046288699e 


- 


01 





3483059E- 


01 





875476530 





959173141 


40 


597846984863 


-iO 


6178727373486 


- 


02 





4986330E- 


01 





885584168 





902043947 


40 


760002136230 


-iO 


5171684897506 


- 


03 





2216994E- 


02 





900596963 





598602624 


39 


372772216797 


-iO 


7821834646176 


- 


02 





3692186E- 


01 





670832741 





857146071 


39 


384689331055 


-iO 


2535248873766 


- 


02 





1560470E- 


01 





644713713 





768504668 


39 


524833679199 


-iO 


4160359967506 


- 


03 





6675268E- 


02 





918124783 





571120689 


40 


427906036377 


-iO 


6540142744786 


- 


01 





2790542E- 


02 





989788008 





989591927 


40 


367130279541 


-iO 


6401911377916 


- 


01 





2049567E- 


01 





949528775 





983221172 


40 


508068084717 


-iO 


8145602047446 




01 





1035770E- 


02 





979746925 





983924594 


40 


537075042725 


-iO 


7174256443986 




01 





6688557E- 


01 





910022901 





939244195 


40 


627620697021 


~iO 


9138841927056 




01 





8095847E- 


02 





858722180 





846152346 


39 


421646118164 


~iO 


8507222682246 




01 





5916540E- 


01 





906494113 





943557886 


39 


419769287109 


-iO 


5661062823496 




04 





10631 11E- 


01 





939688864 





531911062 


39 


733528137207 


-iO 


1750214956706 




01 





3724663E- 


01 





783724495 





988899106 


39 


563480377197 


-iO 


6501793861396 




01 





1242726E- 


02 





990091490 





989958495 


39 


685890197754 


-iO 


7705730199816 




01 





1809454E- 


01 





936621212 





951369186 


39 


650619506836 


-iO 


7902416586886 




01 





1439820E- 


01 





934294618 





935693384 


39 


369716644287 


-iO 


1691374480726 


+ 00 





33488 14E- 


03 





999949714 





999923966 


40 


361137390137 


~iO 


1118765473376 


+ 00 





3443463E- 


02 





998603169 





996169638 


40 


389163970947 


-iO 


6357513484556 




03 





3869480E- 


02 





902893323 





686779751 


40 


273880004883 


-iO 


7099246233706 




01 





1889875E- 


02 





987924322 





995055099 


40 


207023620605 


-iO 


3298169467606 




02 





5157299E- 


01 





870880666 





955577617 


40 


213146209717 


-iO 


2789463102826 




01 





4355001E- 


01 





850281901 





996592453 


39 


904693603516 


-iO 


2241544052966 




01 





2102797E- 


02 





998176457 





998674123 


40 


172363281250 


-iO 


1322104483846 


+ 00 





3454306E- 


03 





999768065 





999390566 


40 


137271881104 


-iO 


8734986931096 




01 





6845250E- 


03 





997220268 





997824221 


40 


025264739990 


-iO 


1622620075946 


+ 00 





1378969E- 


03 





999992333 





999984896 


40 


064556121826 


-iO 


5972184985886 




01 





1994631E- 


03 





998558338 





999194260 


40 


058116912842 


-iO 


7460922352046 




03 





6044292E- 


02 





948958619 





983942653 


39 


997570037842 


-iO 


6594257056716 




01 





1710776E- 


02 





999267906 





999318229 


39 


925022125244 


-iO 


2295515332666 




04 





5157831E- 


04 





999120149 





934566350 



VIII. THE HUSIMI PROJECTION TECHNIQUE FOR OPTICAL DIELECTRIC RESONATORS 

In this section, we will describe the Husimi Projection technique, which allows us to relate a given mode to the 
phase space structures in the SOS. Just as for quantum wavefunctions, for these two-dimensional electromagnetic 
fields we can represent the solutions in real space (the solutions we have been calculating) or, by Fourier transforming 
them, in momentum space. However we are interested in representing the solutions in the phase-space of the problem 
so that we can understand their ray-dynamical meaning, and ultimately in projecting such phase space densities onto 
the SOS which is our standard interpretive tool. Just as in quantum mechanics, we cannot have full information 
about real-space and momentum space at the same time due to the analog of the uncertainty principle, which here is 
reflected by the property of Fourier transforms: 

Ax- Ap> (60) 
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where Ax and Ap are the widths of ^{x) and ^{p) in real and momentum space. Thus our goal is to take the 
solution ipi^) Etnd associate with it a momentum content in some region around each point x, recognizing that 
our resolution in real space is limited by the uncertainty relation. One familiar method for doing this in quantum 
mechanics is the Wigner distribution function|67l|. which preserves exactly certain moments of the wavefunction, but 
has the interpretive problem that it can be negative in some regions of phase space and the practical problem that 
it is typically subject to rapid oscillations. An alternative approach, which is often more useful, was introduced by 
HusimipSlj and can be thought of as a gaussian smoothing of the Wigner distribution. In our context however one 
can think of the Husimi projection as a "windowed" two-dimensional Fourier transform, which is designed to respect 
the uncertainty relation. This is equivalent to projecting on coherent states \z) = \p, x) which are optimally localized 
in both momentum and configuration space, in the sense that the projection saturates the inequality Eq. (|6U|) at its 
lower bound 



Ax = —= = —= An = 
V2k V2 



1 



2fc(To 



(61) 



and equal resolution in both spaces is ensured through the choice of ao. This is a free parameter of the method, and 
for best results it has to be carefully determined based on the domains of variation of the conjugate pair {p,x). The 
real-space representation of a coherent states is given by 



Zxp{x) — ( 2 



1/4 



exp(iA:p • x) exp ( -it\x 



(62) 



where the prefactor ensures normalization. Note that the first exponential factor determines the selection of the 
momentum p (normalized to unity), and the second exponential factor restricts the probe to an isotropic window of 
size Ax = ri/\/2 around x in the configuration space. The Fourier transform is again a localized Gaussian, given by 



1/4 



^xp(p) = — exp(-ifc(p-p) - a;) exp - — \p 



(63) 



Note that both scales Ax and Ap are sharpened as k 
on these states 



oo. We can now construct Husimi distribution by projecting 



d^xz;^jx)^j{x) 



(64) 



Note that this distribution, unlike the Wigner distribution, is positive definite in the phase space {x,p). 

As defined the Husimi distribution is on a four-dimensional phase space of the billiard (restricted to a three- 
dimensional shell because the momentum is normalized to unity) . One can visualize this distribution then as a vector 
field on a grid of size in real space [29^: however based on our earlier discussion we find it more illuminating to 
define a projection of the Husimi distribution onto the surface of section of the billiard. 

For billiard systems, this idea was first carried out in Ref.|l4j. but a different section was used than we arc using 
here. We instead follow an approach similar to the one described in Ref.|23|. The problem is that the boundary of 
the billiard is not a constant coordinate surface of some convenient coordinate system. We therefore use cylindrical 
coordinates to define our Husimi-SOS projection and then use the classical map to map our result onto the boundary. 

The coherent states in cylindrical coordinates take the form: 



where 



1/4 



exp [iprr] exp 



1 

'2^ 



Z^p^i<t>) ^ [zZ^j cxp[ip^{(j)-27rl)]cxp 



(65) 

(66) 
(67) 



The sum on I is necessary to insure periodicity in the </> variable. We define the projection of the full four-dimensional 
Husimi function onto the SOS at constant radius r = Rc hy 



Hi 



lim 
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drZ 



(68) 
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Note that the integration only extends over positive radial momenta to accord with the definition of the SOS which 
only counts trajectories which encounter the boundary in the outgoing direction. We have included in the definition 
Eq. H68|l the limit ry representing the short wavelength limit where k oo. In practice, due to Eqs. I|60|) . I|61(l . 
for any given k value rj is bounded below and is greater than the wavelength A. Violation of this condition would give 
unphysical results. For example, a closed circular billiard would have ip{r = Rc) = and the integrand of Eq. H68() 
would vanish if rj were taken to zero before fc — > oo due to the concentration of the coherent state to a region less than 
a wavelength from the boundary; hence the Husimi-SOS would vanish. Letting r; — > and 1/fc — > while ki] > 1, it 
can be shown that 



lim 



dpr 



d<t) 



drZ 



(r, (j))ip{(j), r) 



where ip^^\Rc, (j)) contains only the wavefunction components with the Hankel's functions of the first kind: 



imcp 



(69) 



(70) 
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The presence of Hankel's functions of only one kind is in accordance with the short wavelength interpretation that 
H^j (nkr) represents incoming and outgoing waves respectively. This interpretation is quickly obscured for components 
m > nkR, which however have vanishing weights in the expansion Eq. H70|) through the physical considerations laid 
out in section (|V1() . 

Using Eq. Ht)9|l and the short wavelength correspondence m ~ nkR^ sin x in Eq- Hfci8|) we obtain 



i?^(0,sinx) 



, X 1/4 oo 



1 — — 00 
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Noting that the integrand in Eq. H71|l is 27r-pcriodic, the integration limits can be extended to infinity and the resulting 
Gaussian integral over can be evaluated analytically, yielding 



iJ^((/),sinx) 



^a™H+(nfci?c)e" 



inkRc (sin x— sin x)4' 



e 2 



(72) 



For optimal resolution in both SOS coordinates, we choose a ~ 1/y/nkRc (Rc can be chosen at any convenient value). 
Finally, in order to calculate the distribution on a section on the boundary r = R{(j)) instead of the circle, we choose 
Rc slightly outside the boundary, and extrapolate the Husimi distribution Eq. H72(l to the boundary using the classical 
equations of motion; i.e. every pair (sinx, 4>) on the circle maps uniquely to a different pair (sinx, 0) on the boundary. 
We can regard this mapping as simply a change of variables so that the Husimi-SOS distribution on the boundary, 
H,p (0, sin x) , satisfies: 



(0((/),sinx),sinx((/',sinx)) = Hip ((/),sinx) 



(73) 



IX. FARFIELD DISTRIBUTIONS 



In typical micro-laser experiments there are two basic data-acquisition modes. In the farfield acquisition mode, the 
CCD camera is used without a lens and aperture as a simple photo-diode, and at each farfield angle 9, the farfield 
emission intensity Ipp{9) is recorded. The farfield emission pattern is one of the few clues we have as to the lasing 
mode of the resonator. Our algorithm is very efficient in calculating all possible farfields achievable with a given 
resonator shape and index of refraction. At a given k which is chosen close to the lasing frequency uo = ck, we solve 
the scattering problem without reference to any quantization condition. The farfields are computed from the external 
wavefunction Eq. H26(l by using the large- argument asymptotic form of the Hankel function 



r 



/(0)« 



(74) 



where we have extracted all the quantities independent of m and (j). The farfield intensity distributions computed 
in this way are very well-behaved and insensitive to the k-valuc used, because the strongly-varying Hankel functions 
drop out of this quantity. 
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FIG. 13: Far-field emission intensity pattern for a "fish" mode, quantized at kR — 19.7392 — 0. 06240* (solid black line) and at 
kR — 21.0210 — 0.06240i (dashed line). The (dash dot) and (dotted) line are the far-field of the unquantized states in between 
with eigenvalues z — —0.4664 + 0.5692j and z = 0.3407 — 0.6472i. Clearly the essential features of the emission patterns are 
calculable with non-quantized modes. In the inset we show the eigenvalues z in the unit circle. 



As seen in Fig. E| the farfields computed are virtually identical to those obtained from the quantized modes as k is 
varied over two level-spacings. Using this simplification the method here can be used to evaluate all possible emission 
patterns for a wide range of dielectric resonators very rapidly. 

Several recent experiments have studied dielectric micro-lasers using an imaging technique for data acauisition |lll 
IsM ls^lsijl : the CCD camera records a magnified image of the intensity profile on the sidcwall viewed from angle 9 in 
the farfield (see Fig. ll4|l . The pixels then can be mapped to sidewall angle (j> via the solution of a simple transcendental 
equation. This yields a two-dimensional plot, called the boundary image-field, where a given data point I{(f>, 9) denotes 
the intensity emitted from sidewall position cf) towards the farfield angle 6. The latter can easily be converted to an 
incidence angle sin x, using SncU's law and basic trigonometry. Hence what is recorded is actually a phase space plot 
of the emitted radiation. 





FIG. 14: a) Experimental setup for measuring simultaneously far-field intensity patterns and boundary images as implemented 
in ref . Isill : the lens and aperture are in the far-field and the schematic is not to scale, b) experimental data on lasing emission 
from the same reference, c) Simulation of the experimental set-up based on calculated modes of the corresponding quadrupolar 
ARC dielectric resonator. 
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FIG. 15: A numerically calculated transverse series of modes based on the bow-tie periodic orbit. Such modes were observed 
in the experiments of ref . |23l . The calculations are performed at nkR ~ 128, n — 3.3, e = 0.16. With these parameters the 
bowtie periodic orbit hits the boundary exactly at the critical angle, a) the fundamental mode, b) the first excited mode c) the 
second excited mode. Note the strong "near-field" fluctuations, particularly in the case of the second order mode in c). In the 
middle row we have their Poincare-Husimi distributions. On the bottom we have for each of the modes a calculated boundary 
image, showing that one can clearly distinguish the three different modes from their boundary-image fields. 



The aperture has an important role of defining a window in the direction space (Asin^), so that a given pixel on 
the camera can be identified up to a diffraction limited resolution with a pair (0,sinx). Mathematically, the effect 
of the lens-aperture combination is equivalent to a windowed Fourier transform of the incident field on the lens|24l|: 
thus it is simply connected to the Husimi distributions we have just discussed [gJ. Note that infinite aperture limit is 
simply a Fourier transform of the incident field and we lose all the information about direction sinx, consistent with 
our intuition with conjugate variables. It has to be emphasized that image data only probes the farfield, and does not 
contain the "near-field" details we would see in a typical numerical solution, nor does it contain information about 
the internally reflected components of the cavity field (see Appendix IBl for further details). On the other hand, it 
contains (with some finite resolution) the same information as the Husimi distribution of the emitting components of 
the field and does allows us therefore to put forward a ray interpretation of the emission and the internal resonance. 

The boundary-image field of a numerically obtained resonance can be calculated as described in Appendix ^ In 
Fig. El we plot the boundary- image fields of three "bow-tie" modes of different transverse excitation number (based 
on a stable periodic orbit with the geometry of a bow-tie, see discussion below). 
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X. QUASI-BOUND MODES AND CLASSICAL PHASE SPACE STRUCTURES 

111 this section, we would like to illustrate the relation between the quasi-bound states and classical phase space 
structures using the results generated by the numerical algorithm described in section ljVII|l and Appendix^ We 
note at the beginning that the real-space plots and the Husimi distributions are constructed using the nonquantized 
scattering eigenvectors, as we described in the previous section. 

Consider first the near-integrable regime in the quadrupolar deformation. In Fig. llBf a). we plot a whispering gallery 
mode of the circle at e = and in Fig. llfif b^ is plotted a state which emits from the highest curvature points (/) = 0, tt. 
Inspection of the respective Husimi distributions shows clearly that the first state is localized on an invariant curve 
sinx = const, and is (nearly) totally internally reflected, because the localization is at sinx > sinxc- Note that the 
second state is a deformed whispering gallery mode, localized on an invariant curve which is no longer a straight line 
in the SOS. 




FIG. 16: Real-space false color plot of regular quasi-periodic solutions at (a) e = 0, and (c) e = 0.03. (b), (d) Husimi projections 
of states (a), (c). The solutions are obtained at nkR — 82 and n = 2. 

Large islands in the SOS (corresponding to stable periodic orbits) support multiple modes, which can be calculated 
approximately by the methods of Gaussian optics [63 • We have already seen a typical sequence of such modes based 
on the bow-tic periodic orbit in Fig. 1151 In general, it's possible to find such series of modes corresponding to any 
stable island of sizable extent in the SOS, with two characteristic spacings, the free spectral range Ak ~ 2tt/L (L is 
the length of the periodic orbit), and the transverse mode spacing, which depends on the eigenvalues of the stability 
matrix |65j- 

The importance of periodic orbits is not limited to stable orbits. There are an infinite number of unstable periodic 
orbits in the SOS. Short periodic orbits, especially those which are least unstable, can make their presence felt in the 
mode structure of the resonators, despite the fact that the methods of Gaussian optics fail for such modesji^- Such 
modes are referred to as "scars" in the quantum chaos literature; they display an enhanced intensity along an unstable 
periodic orbit and have been widely studied [28ll33|. In the numerical studies of Refs.fT^lssj. evidence was found that 
modes can localize not only on the unstable fixed points but on their associated stable and unstable manifolds as well. 
In Fig. [rTT al-ffl we show a series of three states associated with the shortest unstable periodic orbit of the system, 
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FIG. 17: Real-space plots and Husimi distributions of modes which are related to the unstable bouncing ball orbit (scars, see 
text). Superimposed on the SOS are the stable and unstable manifolds of this orbit which seem to play a role in localizing the 
solutions in certain regions of the chaotic component of phase space, (a), (b) represent a simple scar of the two-bounce orbit; 
(c), (d) and (e), (f) are states localized on intersections of the manifolds and appear more chaotic in real space. These solutions 
are found at nkR = 106, e = 0.12 and n = 2.65. 



the unstable bouncing ball orbit. Note that the first mode represents the "fundamental" mode which localizes on the 
fixed point itself (see Fig. I17r b)). The real-space plots Fig. I17f c') and Fig. I17r e) don't show any distinct structure. 
Their respective Husimi plots Fig. I17r d') and Fig. I17r f'l. however reveals that the modes localize on the heteroclinic 
intersections of the stable and unstable manifolds emanating from the unstable bouncing ball fixed points. A similar 
behavior was reported in Ref.|l9l| in the context of a quantum billiard. Finally, at large deformations the spectrum 
contains chaotic modes which cannot easily be associated with any particular classical phase space structure. One 
such mode is plotted in Fig. ^| Note that the support of the mode is entirely in the chaotic portion of the SOS. 
Recalling our arguments in section (|lT|) as to the failure of cikonal methods, it is instructive to note here the complexity 
of the wavefronts and the large portions of the resonator in which no series of parallel wavefronts is discernible. 
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FIG. 18: Real-space lalse-color plot and Husimi distribution of a chaotic mode at a quadrupolar deformation of £ = 0.18 and 
n = 2.65, quantized at kR = 32.6638 - 0.06964i. 

XI. CONCLUSION 

We have considered a general class of optical resonators which are based on deformations of cylindrical dielectric 
resonators. Such resonators are being studied for applications in integrated optics and optoelectronics, and for their 
intrinsic interest as wave-chaotic systems. The basic physics of such resonators can be best understood by application 
of methods from classical and quantum chaos, applicable in the short-wavelength limit. 

We have pointed out the characteristic global breakdown of conventional geometric optics approaches due to the 
transition to chaos in the associated ray dynamics. While real-space ray tracing methods quickly become powerless 
with increasing deformation (degree of chaos), the essential structures are uncovered effectively in the Poincare 
surface of section. Information about important physical properties of deformed dielectric resonators and lasers, such 
as emission characteristics, internal modal distributions, spectra and lifetimes can be extracted from the types of 
ray motion in the equivalent refractive billiard system. However no general analytic technique exists for calculating 
approximately all of the modes of a generic resonator of this type. 

We present an efficient numerical method for the calculation of quasi-bound modes of dielectric cavities. The method 
we are proposing is a hybrid between a poi nt-matching technique |44j and scattering approach to auantization|l6l Il9l| . 
In contrast to existing methods ,.30. l44l l48l| which employ the external scatte ring matrix to extract the quasi-bound 
modes of a dielectric resonator, we consider the internal scattering operator|l5l Il6l Il7j. An important conceptual 
difference with respect to wave-function matching methods is that the internal scattering approach permits the 
identification of a discrete set of internal scattering states at each value of k. This is realized by writing the matching 
conditions in the form of an off-shell eigenvalue problem instead of a linear inhomogcneous equation. Many of the 
physical properties of the modes within a given linewidth of the order of l/nk^R around fco are contained in the 
eigenvectors of the internal S-matrix at fco, and can be extracted without quantizing the mode (i.e. without tuning k 
to satisfy the boundary conditions). The quantized spectrum and the exact quasi-bound modes can be easily accessed 
by an extrapolation technique which requires only two diagonalizations of this S-matrix. In principle variants of this 
technique can be extended to very high wavevectors fc since it scales only the perimeter of the resonator. 

Finally, we present a version of the Husimi projection technique well-suited to dielectric resonators and show that 
the scattering eigenvectors as well as quasi-bound modes are structured by classical phase-space structures. 

In conclusion, the methods and tools presented in this work provide a unified conceptual framework for treating 
dielectric resonators beyond the standard numerical approaches to such electromagnetic problems. 
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APPENDIX A: NUMERICAL IMPLEMENTATION ISSUES 

Although the development in the text using the eigenvalue problem for the truncated S-matrix S{k) is perfectly 
sound, its direct implementation is not very efficient. The numerical problems already surface at the stage of computing 
the truncated scattering matrix S{k) itself. The expression given in Eq. H37|) requires numerical inversion of 5 matrices. 
At a given k, as the number of channels A of the interior matrices designated by index 1 grows beyond A^c, they 
become increasingly more singular. 

As evident from our previous discussion, this ill-conditioning is caused by blindly including evanescent channels 
in the scattering problem. A quick solution is to truncate the matrices at the singularity boundary suggested by 
our ray interpretation and include only A^c = l^ifci^mm] channels of |m|. But this turns out to produce some states 
which do not satisfy the boundary conditions well enough. As noted in|lfi|. some of the evanescent channels Aev 
have to be kept, enough to be able to proceed with our numerical computation and provide the missing (evanescent) 
components of those states which require it. Thus, the properly truncated scattering matrix will have a size of 

The conditioning of S{k) is highly sensitive to the choice of Ae„. Singular value decomposition can be employed to 
calculate the inverses and build up S{k). However, one way to sidestep this issue in favor of a more robust method is 
to trade the eigenvalue problem Eq. for a generalized eigenvalue problem. We rewrite Eqs. (|32I33|) in the form 

A\T) = e''^B|T) (Al) 

where the 2Ntrunc x '2Ntrunc matrices are given by 




A= '"1, 1 ' + , B=\ '-1 " (A2) 



and 



l7) 



(A3) 



This way, the common null-space of the matrices A and B can be removed by existing powerful generalized eigenvalue 
solvers, such as ZGGEV of LAPACK library. As a byproduct we get both the inside and the outside vectors at one shot. 
This method turns out to be more stable than the one based on inversion, and yields good results with arbitrarily 
large Ae^. The numerical problems associated with the regions of evanescent behavior {Rmin < r < Rmax) remain, 

but are tractable for nkR ^ 200 — 300 (this range is larger for modes without evanescent components). 

APPENDIX B: LENS TRANSFORM 

In the experimental imaging system, radiation emanating from the resonator is collected through an aperture and 
after passing through a lens, an image is recorded for a discrete number of angles in the farfield. The resonator is 
placed at the focal plane of the lens, so that the image is effectively formed at infinity. Just in front of the lens the field 
distribution is given by the resonance wavefunction "^(x) which, at the observation (FF) angle 6, can be expressed as: 



*(^) ~ E 7r.H+(fcy^x2 + z2)c™(^+«) (Bl) 

m 

where (p = tan"-'^ ^ ~ F"' ^^'^^ effectively adds a quadratic phase, so that the field immediately behind the lens 
is given by 



^''(a:) = '${x)P{x) exp 



k 2 
'l^r-:X 



2/ 



(B2) 



Here P{x) is the pupil function, which takes care of the effect of the aperture, / is the focal length of the lens and x 
is the position on the lens. The field at the camera is given by propagating this field with the Fresnel propagator [2^. 
which is well-justified as the lens-camera distance is much larger than the wavelength 

*^'(u) = _/ dx^'{x)cxp[i — {u~xf] (B3) 

^Z2 J-oo 22:2 
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FIG. 19: Variables used in the lens transform to generate the boundary image field as measured experimentally. 



Using the expression Eq. (jBl|) for the wavcfunction 

k 

iZ2 



^ig^u) ~ - — e^ ^^2^ ' '2,' 



X r dxP(a;)H+(fcAYz2+x2)e*™(''+^T)e**(^-7) 

J — OO * 



e ^2 



Using the large argument asymptotic expansion of the Besscl functions: 



H+ (fc Jx2 + z\) 



■ exp 



ira I — 

2 4 



(B4) 
(B5) 

(B6) 



fc(x2 + z2)l/2 

Expanding the square root to C'(— )4, i.e. to the same order as the Fresnel approximation and rearranging the terms, 



— e 

iZ2 



-e 



kzi 



1 r,/ \ tm— —ux + — — 4) 

dxP{x)e ^1 ^2 e ^^'1 '2 



(B7) 



(B8) 



The second exponent in the integral is exactly the lens law, so it vanishes. Setting ^ = Ro sin Xm, the intensity 
recorded at the pixel u of the camera at the farfield angle can asymptotically be written as 



/OO 
dxP{x) exp 
'OO 



i — [MRo sinxm — u)x 

Z2 



(B9) 



where M = z^j z\ is the magnification of the lens. For a simple aperture, P(x) is just a rectangle function, so that 
the integral can be performed exactly to yield 



AY,l^H+{kz,)e'"'' sine 



1 , . u . 
— sm 



(BIO) 



where A — and sinc(a;) = sinx/x. Note that in the short-wavelength limit A ^ and ;;^sinc(^) S{x). This 

expression allows us to make predictions based on short-wavelength limit and geometric ray optics, which includes 
effects of diffraction as well. For instance, for a circular cylindrical resonator, the resonances are composed of a 
single angular momentum component m (and its degenerate partner — m). In that case, according to the expression 
Eq. ||BT0)| . 



\^'^{u)f (X |5(sinx,„ - ^y--' ± Sisinx^ + 



(Bll) 
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Note that the imagefield contains only information captured from the farfield distribution. The actual details of the 
resonance in the "nearfield" can be quiet different, due to evanescent contributions close to critical incidence. For 
instance, the points of brightest emission inferred from the imagefield might be shifted due to an "optical mirage" - 
like effect (see Fig. llSfc ')'). The mirage is formed not because of a continuously varying index of refraction but a 
discontinuous interface. 

The imagefield has an interesting connection to the (SOS projected) Husimi distribution. The Husimi distribution 
of the field projected onto the SOS at a distance i? ^ cx) is given by 



^7™H+(fci?)e^'"«e 



im0 -^rf (m-peY 



(B12) 



Comparing with Eq. ljB10|l . we see that the two functions are contain almost the same information. In fact, by 
choosing an aperture which has a gaussian transmittance P{x), one would obtain exactly the same form as Eq. ljB10|l . 
Note that the freedom of smoothing to obtain various phase space distributions which represent the same physical 
system gains here a physical meaning, namely it translates to the choice of optical apparatus (lens, aperture etc.) 
to observe the resonator. This connection was used in reference j6lj to reconstruct from the boundary-image field a 
Husimi function for the emitted radiation which was found to agree well with ray escape simulations. 
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